@vsnever
@mattngc
@CnlPepper
This feature is not prioritized over other issues, but how do you feel about adding this?

I've not had a chance to look at this, but, assuming you haven't, could you make test demo that checks it works correctly with the CSG operations (tests next intersection)?

@CnlPepper
I tried to calculate csg oprations by:

• Torus + Sphere(right) - Sphere(left) - Box(front) (Copper matrial)
• Torus * Cylinder (Glass material)

from matplotlib import pyplot as plt

from raysect.optical import ConstantSF, Point3D, World, d65_white, rotate, translate
from raysect.optical.library.metal import Copper
from raysect.optical.material import Lambert, UniformSurfaceEmitter
from raysect.optical.library import schott
from raysect.optical.observer import PinholeCamera, RGBAdaptiveSampler2D, RGBPipeline2D
from raysect.primitive import (
    Box,
    Cylinder,
    Torus,
    Sphere,
    Union,
    Intersect,
    Subtract,
)


# glass matrial
glass = schott("N-BK7")

world = World()

# Toruses
torus1 = Torus(1.0, 0.5)
torus2 = Torus(1.0, 0.5)

# Spheres
sphere1 = Sphere(0.6, transform=translate(1.0, 0.0, 0.0))
sphere2 = Sphere(0.6, transform=translate(-1.0, 0.0, 0.0))

# Box
sqrt2 = 2 ** 0.5
box = Box(
    Point3D(-1.6, -1.6, -1.6),
    Point3D(1.6, 1.6, 1.6),
    transform=translate(0.0, 2.21, 0.0),
)

# cylinder
cylinder = Cylinder(0.6, 2.0, transform=translate(0.0, 1.0, 0.0))

# Torus1 - Sphere1 + Sphere2 - Box
Subtract(
    Union(Subtract(torus1, sphere1), sphere2),
    box,
    world,
    transform=translate(0.0, 0.0, 0.6),
    material=Copper(),
)

# Torus2 * Cylinder
Intersect(
    torus2,
    cylinder,
    world,
    transform=translate(0.0, 0.3, 0.6),
    material=glass,
)


# floor
Box(
    Point3D(-100, -100, -10),
    Point3D(100, 100, 0),
    parent=world,
    material=Lambert(ConstantSF(1.0)),
)

# emitter
Cylinder(
    3.0,
    100.0,
    parent=world,
    transform=translate(0, 0, 8) * rotate(90, 0, 0) * translate(0, 0, -50),
    material=UniformSurfaceEmitter(d65_white, 1.0),
)

# camera
rgb = RGBPipeline2D(display_unsaturated_fraction=0.995)
sampler = RGBAdaptiveSampler2D(rgb, min_samples=500, fraction=0.1, cutoff=0.01)
camera = PinholeCamera(
    (512, 512),
    parent=world,
    transform=rotate(0, 45, 0) * translate(0, 0, 5) * rotate(0, -180, 0),
    pipelines=[rgb],
    frame_sampler=sampler,
)
camera.spectral_bins = 21
camera.spectral_rays = 1
camera.pixel_samples = 250
camera.ray_max_depth = 10000
camera.ray_extinction_min_depth = 3
camera.ray_extinction_prob = 0.01


# start ray tracing
plt.ion()
for p in range(0, 1000):
    print(f"Rendering pass {p}...")
    camera.observe()
    print()

plt.ioff()
rgb.display()
plt.show()

I think the implementation of the next_intersection() method is inconsistent with its intended behaviour. The description of this method says:

    cpdef Intersection next_intersection(self):

        if not self._further_intersection:
            return None

        # this is the 2nd intersection
        self._further_intersection = False
        return self._generate_intersection(self._cached_ray, self._cached_origin, self._cached_direction, self._next_t)

The method is used only in the CSGPrimitive._identify_intersection() method.

As a result, extra intersections are generated in some cases. The output of this script:

from raysect.optical import Point3D, Vector3D, World, translate, InterpolatedSF
from raysect.optical.material import Dielectric
from raysect.primitive import Union, Torus, Cylinder
from raysect.optical.loggingray import LoggingRay


glass = Dielectric(index=InterpolatedSF([10, 10000], [1., 1.]),
                   transmission=InterpolatedSF([10, 10000], [1., 1.]),
                   transmission_only=True)

world = World()

torus = Torus(1.0, 0.5)

cylinder = Cylinder(1.0, 1.0, transform=translate(0, 0, -0.5))

union = Union(torus, cylinder, material=glass, parent=world)

ray = LoggingRay(origin=Point3D(2., 0, 0), direction=Vector3D(-1, 0, 0))
ray.trace(world)
for intersection in ray.log:
    print(intersection.hit_point, intersection.exiting)
    print()

is

Point3D(1.5, 0.0, 0.0) False

Point3D(1.4999997776219791, 0.0, 0.0) False

Point3D(-1.0, 0.0, 0.0) True

Point3D(-1.5, 0.0, 0.0) True

The intersection at Point3D(-1.0, 0.0, 0.0) is an extra one because the ray is not exciting the primitive at (-1, 0, 0).

Oh this looks fun! I'd need to check this over with a fine tooth comb.... has the comment from @vsnever regarding the next_intersection() implementation been addressed?

@CnlPepper Thank you for commenting on this stale PR.
Since the implementation of the 4-times intersection is very specific to the torus primitive, I wonder whether we should modify other related codes.
Also, I haven't yet found a workaround for the issues mentioned above.
Could you please help make some changes to resolve them?

The next_intersection() implementation just needs to return each interection found in order. For most of our current primitives (except mesh) there are only two intersections. So hit() returns the first, with the second cached. A subsequent call to next_intersection() returns the second, cached one. For the torus you'll need to cache up to 3 addition intersections on the object and return them one by one, in ray distance order for each call of next intersection.

Thank you for giving me informative advice!
I will try to implement the up to 3 hit-point caching functionality, and confirm if the above issue will be resolved.

The result of executing the above script that @vsnever proposed is as follows:

Point3D(1.5, 0.0, 0.0) False

Point3D(1.4999997776219791, 0.0, 0.0) False

Point3D(-1.5, 0.0, 0.0) True

It seems to be an appropriate behavior.
I am wondering a bit why the second one appeared, and the hit points were not symmetrically represented.

Note that some lines are still incompatible with Cython 3.1, and we should address them at #447 if this PR is merged beforehand.

I would like to modify the document about the primitive page by rendering the top image, adding a torus primitive.
Does anyone know where the script that renders its image is located? Or should I write it from scratch?

I would like to modify the document about the primitive page by rendering the top image, adding a torus primitive. Does anyone know where the script that renders its image is located? Or should I write it from scratch?

I think that was a one-off script I created in the early Raysect development. I have a quick look but couldn't find it. If you recreate it, feel free to add it to demos/primitives if you like.

I'm pretty sure the glass material I used was N-Bk7. You'll need to hover the primitives slightly off the ground plane (1e-6 m or so) to prevent surface fighting with the ground plane.

I created a primitive image below. I will replace the previous image and add the script file, rendering it.

When executing the generate_meson_files.py script, almost all meson.build files were modified, which seems to reorder the elements of the source files lists.
The pathlib.Path.iterdir() here:

Interesting, I'd not encountered that thus far. I'll make all the lists sorted to avoid the issue in the future. Thanks for spotting that.

I believe this PR is ready to merge!